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We present results of 3D numerical simulations using a finite difference code featuring fixed mesh 
refinement (FMR), in which a subset of the computational domain is refined in space and time. We 
apply this code to a series of test cases including a robust stability test, a nonlinear gauge wave and 
an excised Schwarzschild black hole in an evolving gauge. We find that the mesh refinement results 
are comparable in accuracy, stability and convergence to unigrid simulations with the same effective 
resolution. At the same time, the use of FMR reduces the computational resources needed to obtain 
a given accuracy. Particular care must be taken at the interfaces between coarse and fine grids to 
. . . avoid a loss of convergence at higher resolutions, and we introduce the use of "buffer zones" as one 

' resolution of this issue. We also introduce a new method for initial data generation, which enables 

' higher-order interpolation in time even from the initial time slice. This FMR system, "Carpet", is 

, a driver module in the freely available Cactus computational infrastructure, and is able to endow 

(N . generic existing Cactus simulation modules ("thorns") with FMR with little or no extra effort. 

PACS numbers: 04.25.Dm 04.20.-q 04.70.Bw 

cni : 

■ I. INTRODUCTION 

> ; 

^3 ' Currently many researchers are involved in efforts to predict the gravitational waveforms emitted by the inspiral 
' and merger of compact binaries. The direct numerical simulation of these binary sources has been regarded as 
the best, if not only, way to obtain information about the waves emitted during the merger event itself. Many 
T— I numerical simulation methods are currently in use, and among the most popular of these is the use of finite difference 
approximations to Einstein's equations on a uniform mesh or grid. 

It is often the case in simulations of physical systems that the most interesting phenomena may occur in only a 
O subset of the computational domain. In the other regions of the domain it may be possible to use a less accurate 
approximation, thereby reducing the computational resources required, and still obtain results which are essentially 
^ similar to those obtained if no such reduction is made. In particular, we may consider using a computational mesh 
which is non-uniform in space and time, using a finer mesh resolution in the "interesting" regions where we expect 
it to be necessary, and using a coarser resolution in other areas. This is what we mean by mesh refinement (MR). In 
' ^ the case of large-scale simulations, where memory and computing time are at a premium, the effective use of MR 
can allow for simulations achieving effective (fine-grid) resolutions which would be, on a practical level, impossible 
otherwise. This also implies that some systems requiring a large dynamic range of resolution may only be possible 
■ - - to simulate using mesh refinement. One hopes that the use of fewer total grid points in an efficient mesh refinement 
application, as opposed to a uniform, single grid or unigrid application, may result in a saving of computing time 
required to arrive at a solution. It is always the case, however, that the reduction of memory requirements afforded 
by mesh refinement, especially in 3D, could allow for highly accurate simulations which simply may not fit in the 
memory of a given computer in a unigrid setting. 

In many simulations it may be desirable for the computational mesh to adaptively refine in response to some 
criterion such as the local truncation error. Such adaptive mesh refinement (AMR) systems provide, ideally, a very high 
degree of computational efficiency for a given problem. AMR schemes have been in use for decades 1 1], particularly 
in the fluid d^jmamics community jUtUb- Applications to to astrophysical systems have been underway for quite some 
time (e.g., f?, 5,"^, "7]). In numerical relativity, Choptuik f^, 9] introduced many researchers to AMR through his results 
using a ID code for studying critical phenomena. Since the early 90's efforts have been underway in the numerical 
relativity community to develop and e mploy A MR applications in two and three dimensions, in studies of waves 
ESIllllllllllSliSl], critical collapse l^llllllliSi], the initial data problem inhomogeneous cosmologies 
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Schwarzschild black holes characteristic methods fl^, and a binary black hole system fl^. In 

several of these cases, the use of mesh refinement provided not only a more efficient means of computing a solution, 
but even allowed the authors to obtain solutions which would not have been possible via standard unigrid methods, 
given existing computing technology at the time. 

One step in the development of an AMR code is the ability to handle regions of different resolution which are 
known in advance — fixed mesh refinement (FMR). Thus we view our present work as in some sense a precursor 
to full 3D AMR simulations. However, the FMR code may be sufficient in its own right, as full AMR may not be 
necessary for problems where the region of interest is well known beforehand. 

An important point should be made regarding the benefits of mesh refinement for reducing storage requirements 
(and computational costs). It is well known that pseudospectral collocation methods offer exponential convergence 
as the number of collocation points is increased (for sufficiently smooth functions, which are present in many systems 
of interest to numerical relativists), and such methods are able to yield extremely high accuracies via simulations 
which require storage afforded by a typical desktop computer |26, 27, .28.. .29,, 30]. Since FMR makes use of finite 
difference methods, we can obtain at best only polynomial convergence as we increase the resolution of all grids. 
Our choice of finite difference methods and FMR are based on a choice of physical systems of interest, and on an 
observation. 

One reason we are motivated to use mesh refinement is that we are interested in systems which are "non-smooth" 
by the standards of pseudospectral methods: hydrodjmamics, and gravitational collapse. In hydrodynamics, the 
formation of discontinuities — shocks — from smooth initial data is a generic phenomena. In gravitational collapse, 
features may appear on smaller and smaller spatial scales, requiring a means of resolving these features as they 
appear, either via a truly adaptive algorithm (in which, e.g., the truncation error is used as a refinement criterion 
111]) or a "progressive" mesh refinement system in which nested grids are "turned on" at smaller radii and higher 
resolutions as the simulation proceeds. 

The observation mentioned above is that unigrid finite difference codes are generally regarded as being signifi- 
cantly easier to develop than pseudospectral codes. Our interest, from a code-development standpoint, has been to 
provide an infrastructure whereby existing unigrid codes can be "endowed" with mesh refinement in a way which is 
somewhat automatic, such that the developer of the original unigrid code is spared the details of implementing mesh 
refinement. If the introduction of mesh refinement does not significantly alter the dynamics of the system, then one should 
be able to obtain results comparable to a high resolution unigrid nm via an FMR nm with appropriately placed fine 
grids which share the resolution of the imigrid run. In fact, this is the hope of this paper, and the criterion by which 
we evaluate our results obtained by FMR: ideally, the use of mesh refinement should produce results of comparable quality 
to a corresponding unigrid run, in terms of stability, accuracy and convergence behaviour. Thus our mesh refinement 
infrastructure could, in applicable cases, provide a service to the community of researchers who commonly develop 
unigrid finite difference codes, by providing these researchers a means by which to achieve more accurate results 
than their current computer allocations allow. It is in this spirit that we are making the FMR system, called Carpet 
(authored by Erik Schnetter, with refinements offered by several others), freely available as a driver thorn of the 
open-source Cactus computational infrastructure l3lll3Z,l33,l33.l3al . 

II. OVERVIEW OF OUR FMR METHOD 

A. Cactus and mesh refinement 

Cactus is an application framework that provides some computational infrastructure such as parallelisation and 
disk I/O, and lets users specify the modules, called "thorns", which contain all the "physics". The main difference 
between an application framework and a library of subroutines is control inversion, meaning that it is Cactus calling 
the users' routines, while the main program is part of Cactus. The part of Cactus that controls the creation of initial 
data and the time stepping is called the driver, which is a thorn that interacts with the Cactus scheduler in order to 
determine which routines are applied to which grid at what time. Control inversion has the important advantage 
that, by replacing the driver, one can change a Cactus application from unigrid to mesh refinement without rewriting 
any of the users' thorns. 

In practice, the ability to use the same code with both imigrid and mesh refinement drivers places restrictions on 
the implementation of a routine. For most evolution thorns these modifications are at most minor. Some thorns, 
particularly those solving elliptic equations, require more substantial alterations. In most circumstances the restric- 
tions imposed are only technical points, and it should be simple for any new code to be implemented to work with a 
unigrid or mesh refinement driver interchangeably. Our experience is also that analysis tools for e.g. wave extraction 
or apparent horizon finding continue to work almost without changes. 
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B. Mesh refinement method 

The mesh refinement driver that we use is called Carpet and is available together with the application framework 
Cactus. It uses the Berger-Oliger approach |1], where the computational domain as well as all refined subdomains 
consist of a set of rectangular grids. In particular, we base our scheme on the so-called "minimal Berger-Oliger" 
setup popularised by Choptuik |36]. In this simplified version of Berger-Oliger, the grid points are located on a 
grid with Cartesian topology, and the grid boundaries are aligned with the grid lines. In our version, we also allow 
fine grid boundaries to occur in between coarse grid points, as shown in Figure Q] (Note that this definition still 
allows, e.g., spherical coordinate systems.) Furthermore, there is a constant refinement ratio between refinement 
levels (described below). 

We use the following notation. The grids are grouped into refinement levels (or simply "levels") L*^, each containing 
an arbitrary number of grids G'^y. Each grid on refinement level k has the grid spacing (in one dimension) Ax*^. The 
grid spacings are related by the relation Ax^ = Ax'^^^/Nrefme with the integer refinement factor Nrgfine- '^ri example is 
shown in Figure In what follows we will assume that Nj-ggng is always set to 2. The base level I? covers the entire 
domain (typically with a single grid) using a coarse grid spacing. The base level need neither be rectangular nor 
connected. The refined grids have to be properly nested. That is, any grid G'^j must be completely contained within 

the set of grids L^~^ of the next coarser level, except possibly at the outer boundaries. 

Although Cactus does allow both vertex and cell centred grids, current relativity thorns only use vertex centring. 
Hence Carpet currently only supports vertex centred refinement, i.e. coarse grid points coincide with fine grid points 
(where fine grid points exist). 
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FIG. 1: Base level Gq and two refined levels G| and G?, showing the grid alignments and demonstrating proper nesting. 

The times and places where refined grids are created and removed are decided by some refinement criterion. The 
simplest criterion, which is also indispensable for testing, is manually specifying the locations of the refined grids 
at fixed locations in space at all times. This is called fixed mesh refinement. A bit more involved is keeping the same 
refinement hierarchy, but moving the finer grids according to some knowledge about the simulated system, tracking 
some feature such as a black hole or a neutron star. This might be called "moving fixed mesh refinement". Clearly 
the most desirable strategy is an automatic criterion that estimates the truncation error, and places the refined grids 
only when and where necessary. This is what is commonly imderstood by adaptive mesh refinement. Carpet supports 
all of the above in principle, but we will only use fixed mesh refinement in the following. It should be noted that 
automatic grid placement is a non-trivial problem (see, e.g., I37..38.I'). 

C. Time evolution scheme 

The time evolution scheme follows that of the Berger and Oliger fl'l AMR scheme, in which one evolves coarse 
grid data forward in time before evolving any data on the finer grids. These evolved coarse grid data can then 
be used to provide (Dirichlet) boundary conditions for the evolution of data on the finer grids via prolongation, i.e. 
interpolation in time and space. This is illustrated in Figure|2 For hyperbolic systems, where a Courant-like criterion 
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holds, a refinement by a factor of Nrefine requires time step sizes that are smaller by a factor Nj-ggne/ and hence Nrefine 
time steps on level k + 1 are necessary for each time step on level k. At time steps in which the coarse and fine 
grids are both defined, the fine grid data are restricted onto the coarse grid (via a simple copy operation) after it has 
been evolved forward in time. If there are more than two grid levels, then one proceeds recursively from coarsest 
to finest, evolving data on the coarsest grid first, interpolating this data in time and space along boundaries of finer 
grids, evolving the finer grid data, and restricting evolved data from finer to coarser grids whenever possible. This 
is illustrated in Figure|3] 



• — • — • — • 



/M , , , , 



FIG. 2: Schematic for the prolongation scheme, in 1 + 1 dimensions, for a two-grid hierarchy. The large filled (red) circles represent 
data on the coarse grid, and smaller filled (green) circles represent data on the fine grid. The arrows indicate interpolation of coarse 
grid data in space and time, necessary for the boundary conditions on the fine grid (explained in section ing . 
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FIG. 3: Schematic for the time evolution scheme, in 1 + 1 dimensions, for a two-grid hierarchy. The large filled (red) circles 
represent data on the coarse grid, and smaller filled (green) circles represent data on the fine grid. The algorithm uses the 
following order. 1: Coarse grid time step, 2 and 3: fine grid time steps, 4: restriction from fine grid to coarse grid. Since the fine 
grid is always nested inside a coarse grid, there are also coarse grid points (not shown) spanning the fine grid region (at times 
when the coarse grid is defined) at the locations of "every other " fine grid point; the data at these coarse grid points are restricted 
(copied directly) from the fine grid data. 



For time evolution schemes that consist only of a single iteration (or step), the fine grid boundary condition needs 
to be applied only once. Most higher-order time integrations schemes, such as Runge-Kutta or iterative Crank- 
Nicholson, are actually multi-step schemes and correspondingly require the fine grid boundary condition to be 
applied multiple times. If this is not done in a consistent manner at each iteration, then the coarse and the fine grid 
time evolution will not couple correctly, and this can introduce a significant error. We explain this in more detail in 
AppendixlAl 

There are several ways to guarantee consistent boimdary conditions on fine grids. Our method involves not pro- 
viding any boundary condition to the individual integration substeps, but instead using a larger fine grid boundary, 
as demonstrated in Figure|31 That is, each of the integration substeps is formally applied to a progressively smaller 
domain, and the prolongation operation re-enlarges the domain back to its original size. Note that this "buffering" 
is done only for prolongation boundaries; outer boundaries are handled in the conventional way. Also, this is done 
only for the substeps due to the time integration scheme, so that the prolongation is applied at fine grid times when 
there is no corresponding coarse grid time. Note also that the use of buffer zones is potentially more computationally 
efficient. 

We emphasise that the use of these buffer zones is not always necessary. To our knowledge the buffer zones are 
necessary only when the system of equations contains second spatial derivatives, and a multi-step numerical method 
is used for time integration. This issue arises for the BSSN system discussed below. We also give a simple example 
using the scalar wave equation in section HV Al and appendix IaI 



D. Inter-grid transport operators 



As described above, the interaction between the individual refinement levels happens via prolongation and restric- 
tion. For prolongation. Carpet currently supports polynomial interpolation, up to quadratic interpolation in time, 
which requires keeping at least two previous time levels of data. It also supports up to quintic interpolation in space, 
which requires using at least three ghost zones. We usually use cubic interpolation in space, which requires only two 
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after prolongation 
final substep 
second substep 
first substep 
previous time 

FIG. 4: Schematic for the "buffering" during time integration. Shown is the left edge of a refined region, which extends further 
to the right, which is integrated in time with a 3-step ICN method. At the filled points (in the interior), time integration proceeds 
as usual. The empty points (near the boundary) are left out, because no boundary condition is given during time integration. A 
prolongation after the time integration fills the empty points again. This whole scheme corresponds to either of the steps labelled 
2 and 3 in Figure|3| 

ghost zones. (Quadratic interpolation in space introduces an anisotropy on a vertex centred grid.) For restricting. 
Carpet currently uses sampling (i.e., a simple copy operation). These transport operators are not conservative. Since 
our formulation of Einstein's equation (see below) is not in a conservative form, any use of conservative inter-grid 
operations offers no benefit. However, the transport operators can easily be changed. (For more discussion of the 
situations where conservative inter-grid operators are useful or not, see L39ll .) 

E. Initial data generation 

Initial data generation and time evolution are controlled by the driver Carpet. Initial data are created recursively, 
starting on the coarsest level L^. This happens as follows: On refinement level L^, the initial data routines are called. 
This fills the grids on this level. Then the refinement criterion is evaluated (which might be nothing more than a fixed 
mesh refinement specification). If necessary, grids on a finer level L^^^ are created, and initial data are created there, 
and on all finer levels recursively. Then, the data from level L^'^^ are restricted to level to ensure consistency. 

In many cases, the initial data specification is only valid for a single time t = 0, such as when using a time- 
symmetric approach, or when solving an elliptic equation. However, for the time interpolation necessary during 
prolongation (see above), it may be necessary to have data on several time levels. One solution is to use only lower 
order interpolation during the first few time steps. We decided instead, according to the Cactus philosophy, that the 
data that are produced during the initial data creation should in principle be indistinguishable from data produced 
by a time evolution. Hence we offer the option to evolve coarse grid data backwards in time in order to provide 
sufficient time levels for higher order interpolation in time at fine grid boundaries. This ensures that no special case 
code is required for the first steps of the time evolution. 

This initial data generation proceeds in two stages. First the data are evolved both forwards and backwards in time 
one step, leading to the "hourglass" structure illustrated Figure|5| This evolution proceeds recursively from coarsest 
to finest, so that all data necessary for time interpolation are present. Note that this would not be the case if we 
evolved two steps backwards in time, as there would not be enough data for the time interpolation for the restriction 
operation between these two steps. 

In the end we must provide initial data only at times preceding the initial time t = Q; i.e., the hourglass structure of 
Figure|5]is invalid as an initial data specification in Cactus. Therefore we perform in the second stage of this scheme 
one additional step backwards in time on each level, leading to initial data at the times t^, to — At^, and t^ — 1/S.t^ on 
each level L^. 

III. PHYSICAL SYSTEM 

The set of equations we solve are described in detail in |40], and although we briefly review the material here, 
we suggest interested readers refer to the prior publication. The evolution system is that of Shibata-Nakamura |41j 
and Baumgarte-Shapiro li^l , the so-called BSSN formulation. The physical quantities present in a typical ADM |43] 
evolution are the 3-metric y,^ and the extrinsic curvature Kjj. In the BSSN formulation, one instead evolves a different 
set of variables: Kjj is decomposed into its trace K and its trace-free part 




(1) 
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t = +AtO 

t = H-At"/ 2 

t = +AtO/4 

- - - - - - t = 

t = -AtO/4 
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FIG. 5: Schematic for initial data scheme, in 1 + 1 dimensions. Our use of quadratic interpolation in time requires three time 
levels of coarse grid data in order to provide boundary data for evolution on fine grids. To achieve this from the beginning of 
the evolution (without the use of a known continuum solution with which to "pre-load" these levels), we evolve our initial data 
(defined at t = 0) both forwards and backwards one step in time. In this way, three time levels of coarse grid data are always 
available to provide boimdary data along the edges of fine grids. The data at various times are denoted by fractions of the time 
step Af*^ on the base grid. The coarsest grid is shown as by a solid (red) line, a finer grid by a long-dashed (green) line, and a still 
finer grid by a dotted (blue) line. (We perform some additional backwards evolution as well, which we describe in the main text. 
The essence of the scheme, however, is given here.) 

and one applies a conformal transformation, 

fij = e-'*r,r (2) 

We choose cp such that the determinant of fij, denoted by y, is 1. 

Thus our evolved quantities are related to the ADM and physical quantities by 

4^ = ^log(r) (3) 

fij = (4) 
K = yiKij (5) 

Aij = e-^t'iKij - ^YijK) (6) 
One also creates a new evolved variable P', defined as 

r' := r]^iK (7) 

The evolution equations for these variables are given by 

^t(^) = -laK + l3%cp + ld,,l3'' (8) 
6 6 

dtfij = -2aAij + I3%fij + fikdjl3'' + fjkd,l3'' -^f.jdkli'' (9) 

dtK = -D'Dia + a{AijA'i + ^K^) + 13'diK (10) 
dtAij = e-^'t'[-DiDja + aRjjf^ + a{KAij - 2AikA)) 

+A^jdil3^ + AkidjH'' - ^A,jdkl3'' (11) 

, Al<i 
'jk 



-IdjaA'j + 2a{r'^A''' - "^f'djK + 6A'idjcp) 



-djil3%f'' - f^'dkli' - ^'dkl3' + If'dkP"), (12) 

where D, is the covariant derivative corresponding to the 3-metric y/y, Rjj is the three dimensional Ricci tensor, and 
the "TF" superscript denotes the trace-free part of the enclosed expression. 

The gauge conditions are given as follows. The lapse a is chosen as one of the Bona-Masso jiSl family of slicings, 

dta= -a^f{a){K-Ko) (13) 
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where Kq = K{t = 0), and f{a) is chosen to give us either harmonic slicing (/(a) = 1) or a "1 + log" slicing 
(/(a) = 2/a). In this paper the shift will be held constant at the analytic value in all cases. 

We use the code described in [45, 46] and refer there for details of the finite differencing scheme used. 

For certain problems a small amount of artificial dissipation is useful. We use a dissipation of the form 

dtU = -eh"^ {dy:xxx + dyyyy + dzzzz)u (14) 

with the grid spacing h, which is described in l4^ . Although this dissipation operator is typically employed only in 
systems with first order derivatives in time and space, we find that its use in the BSSN system (which has first order 
derivatives in time, and second order derivatives in space) is effective at reducing high-frequency oscillations (i.e., 
noise) in the simulations, but has little effect on the overall convergence behaviour. 

IV. TESTS 

For simplicity we will present tests using only two levels containing one grid each, which we will refer to as the 
"coarse grid" and the "fine grid". The fine grid is a box contained in the larger coarse grid box, with the fine grid 
having a mesh spacing (in space and time) of half that of the coarse grid. The only limitation on the number of 
grids is the available computational resources, and we have successfully performed tests with up to 24 levels of 
refinement. 

One of the principal criteria we use to evaluate the effectiveness of the FMR scheme is the requirement of second- 
order convergence in the limit as the mesh spacing goes to zero. Thus we run a given simulation many times at 
different resolutions. In our examples we compare against the exact solution and check the convergence of the 
solution error. We show such tests only for the data on the coarsest grid, because the restriction operator ensures 
that the coarse grid and fine grid data are identical. 

A. Wave Equation: Periodic Boundaries 

We tested our code with a simple wave equation rn flat space in Cartesian coordinates using several different 
kinds of initial data and boundary conditions. The first such test was that of sinusoidal plane waves in a 3D box with 
periodic boundary conditions. From a code-development standpoint, we simply took an existing set of subroutines 
for solving the wave equation in parallel and ran using the FMR driver instead of the usual unigrid driver. The 
formulation of the wave equation we used was a single equation with second order derivatives in both time and 
space, i.e. 

dttU = d'dfU, (15) 

which we solved using a leapfrog-like scheme. Second order convergence was found, as shown in Figure|6l 

An alternative formulation of the wave equation uses second order derivatives in space but only first order rn 
time. We write it rn the form 

dtu = V (16) 
dfV = d'diU (17) 

This formulation is comparable to the ADM (and BSSN) formulations of the Einstein equations. When this formula- 
tion of the scalar wave equation is evolved using ICN integration without buffer zones, the result is only first-order 
convergent, as shown using a one dimensional example in Figure|7| The same formulation evolved with buffer zones 
converges even at extremely high resolution as shown in FigurelS] 

Having demonstrated the need for careful handling of refinement boimdaries, and having introduced buffer zones 
as an effective approach to that, we use buffer zones rn all the remaining tests discussed in this paper. 

B. Wave Equation: Gaussian pulse 

We consider a Gaussian pulse that crosses a mesh refinement boundary, travelling from the fine into the coarse 
region. This is supposed to mimic the case of gravitational waves propagating from fine, inner grids and radiating 
out into coarser grids. 
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FIG. 6: Result at 50 crossing times in the evolution of the wave equation with periodic boundaries, for two-level runs at three 
different resolutions (1/40, 1/80, and 1/160). Here we show only data for the coarse grid; the finer grid covers the box x,- 6 
[—0.25,-1-0.25]. The left panel shows slices along the x-axis of the (3D) function u found numerically, as well as the analytic 
(i.e. continuum) solution. (In this graph the highest resolution data and the analytic solution appear to lie on top of each other) 
One can see a resolution-dependent phase shift in u. The right panel shows the solution error, defined as the difference between 
the u obtained numerically and the continuum solution. These error graphs have been scaled to demonstrate the second order 
convergence of the numerical results. 
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FIG. 7: Scaled solution errors for the 1-D scalar wave equation at 0.75 crossing times, discretised ivithout buffer zones. The left 
panel shows the error in the solution u, the right panel in its time derivative v. At low resolutions, i.e. from 1/20 up to about 
1/1280, the scheme seems to be second-order convergent. However, at higher resolutions it becomes clear that this is not the 
case. At lower resolutions, the refinement boundaries are visible as small discontinuities in the error of o at x = ±0.25. At higher 
resolutions, the discontinuity develops an oscillating tail and propagates through the simulation domain. It is instructive to see 
that 8 convergence test levels were necessary to see this behaviour numerically. Compare this against Figure|S| 



We use an effectively one dimensional domain (planar symmetry in 3D) with x E [—0.5, -1-0.5], and a coarse grid 
resolution of Ax^ = 1/100. The region x > is refined by a factor of 2. The Gaussian pulse starts in the refined 
region and travels to the left. Figure |5] shows the pulse after it has crossed the interface, and compares the result to 
two unigrid simulations. The errors that are introduced at the refinement boundary are very small and converge to 
second order. In particular, at this resolution about 10^^ of the original pulse is reflected by the refinement boundary. 



C. Wave Equation: 1/r, with excision 



Another test we performed was the solution of the wave equation <15t using initial data which scaled inversely 
with r, the distance from the origin. This was conceived as a useful step before moving on to black holes |48], as the 
"puncture data" |49] we use for the black holes has elements which scale as 1/r in the vicinity of the puncture. Other 
systems with isolated central masses may also be expected to have elements which scale in a similar fashion. The 1 /r 
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FIG. 8: Scaled solution errors for the 1-D scalar wave equation at 0.75 crossing times, discretised ivith buffer zones. The left 
panel shows the error in the solution u, the right panel in its time derivative v. The refinement boundaries are visible as small 
discontinuities in the error of v at x = ±0.25. The results are second-order convergent up to extremely high resolutions. It can 
also be seen that the resolution 1 /20 is clearly not yet in the convergence regime. 




FIG. 9: The left panel shows a Gaussian pulse that has travelled through a refinement boimdary located at x = 0. The region 
X > is refined. The right panel shows the difference to the analytic solution, and shows two unigrid runs for comparison. The 
resolutions of the unigrid runs correspond to the coarse and the fine regions in the mesh refinement run. As expected, the error 
in the run with refinement is in between the errors of the two unigrid runs. The reflected part of the pulse is very small; it has an 
amplitude of only 10^^ and is thus much smaller than the discretisation error in the pulse itself. 



data are a static solution of the wave equation in 3D, and are compatible with the standard Sommerfeld outgoing 
boundary condition. To handle the singularity at the centre, we "excise" the centre of the computational domain by 
choosing some inner boundary at a finite excision radius from the centre and filling the interior region with prescribed 
data. Thus in this test we also have a test of the use of FMR in the presence of excision, which is commonly used for 
the interiors of black holes in analogous simulations. We perform no evolution within this excised region. 

We use a full 3D grid with G [—1, +1] and coarse grid resolutions of Ax^ = 1/32 and Ax^ = 1/64. The region 
X, £ [—0.5, +0.5] is refined by a factor 2, and the region | < 0.125 is excised. Graphs of the error at two different 
times in the evolution are shown in FigurelTUI which also shows corresponding unigrid runs for comparison. We see 
that the solution is fully convergent, and similar to the corresponding unigrid results in the region of refinement. It 
is interesting to note that even outside this region, the FMR and unigrid results are very similar for the "transient" 
shown in the left frame, however the late time ("stationary") behaviour shown in the right frame reveals a notable 
difference between the FMR and unigrid results outside the refined grid. 

Having demonstrated the existence of convergent solutions of the wave equation for oblique angles of incidence 
to refinement boundaries (in the 1/r case), and convergent solutions in which the amplitude of the reflected wave 
is significantly (roughly three orders of magnitude) less than the amplitude of the transmitted wave, we now move 
on to full solutions of Einstein's equations. For more detailed calculations of reflection and transmission effects at 
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FIG. 10: Solution error for evolutions of 1/r data at times t = 0.625 (left frame) and t = 5.0 (right frame), using radiation 
(Sommerfeld) boundary conditions. Here we show only the base grid; the refined grid covers the domain x,- 6 [—0.5, +0.5]. The 
region | < 0.125 is excised. We see in the left frame that an initial transient moving outward from the refined region is still 
resolved in a second-order-convergent manner after it passes through the refinement boundary. At later times, the error settles 
down to the profile shown in the right frame. As expected, the FMR results compare favourably with the unigrid results in the 
region covered by the refined grid. (Note that, outside the refined region, the errors shown for the FMR results have effectively 
been multiplied by a factor of 4 relative to the comparable unigrid results, giving the impression of greater disagreement than is 
really present.) 

refinement boundaries see 



D. Robust Stability 

We have applied a robust stability test jsollsills^l to the BSSN formulation. This is a test that is meant to supple- 
ment and extend an analytic stability analysis, especially in cases where such an analysis is difficult or impossible 
(because the equations contain too many terms). A numerical test has the advantage that it tests the complete 
combination of evolution equations, gauge conditions, boundary conditions, as well as the discretisation and the 
implementation, and (in our case) the mesh refinement scheme. Thus, while a numerical test is not as reliable as an 
analytically obtained statement, it is able to cover more general cases. 

The robust stability test proceeds in three stages of increasing difficulty: 

Stage I: Minkowski (or some other simple) initial data with small random perturbations. The simulation domain is 
a three dimensional box with periodic boundary conditions. The perturbations should be linear, so we chose a 
maximum amplitude of 10^^*^. The periodicity means that there is effectively no boundary, so that this stage is 
independent of the boundary condition. A code is deemed to be robustly stable if it shows at most polynomial 
growth and if the growth rate is independent of the grid resolution. This is different from other definitions of 
stability, where exponential growth is often deemed to be stable if the rate of exponential growth is independent 
of the grid resolution. 

Stage II: The same as stage I, except that the boundaries in the x direction are now Dirichlet boundaries. In addition 
to the noise in the initial data, noise is also applied to these Dirichlet boundaries. This tests the consistency of 
the formulation with the boundary conditions, but without the complications of edges and corners. 

Stage III: The same as stage II, except that there are now Dirichlet boundary conditions in all directions. This tests 
whether edges and corners are handled correctly by the combination of the formulation and boundary condi- 
tions. 



In accordance with 15211 we chose an effectively one dimensional, planar symmetric domain that extends almost 
only in the x-direction with x e [—0.5, +0.5]. The domain has only 3 grid points in the y and z directions. Thus we 
have to omit stage III of the test here. We used a resolution of 1 /50, and refined the centre of the domain. Although 
the domain was thus essentially one dimensional, the simulation was performed with the full three dimensional 
code. 

Figure^Jshows the Loo-norm of the Hamiltonian constraint versus time for 1000 crossing times for the stages I 
and II for three different resolutions. All test were run without artificial dissipation. The results show that our 
implementation is robustly stable. 
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FIG. 11: Loo -norm of the Hamiltonian constraint vs. time for the robust stability test. The left hand side graph shows stage 1 where 
all boimdaries are periodic. The right hand side graph shows stage II, which has noisy Dirichlet boundaries in the x direction. 
The domain is essentially one dimensional with x 6 [—0.5, +0.5], and periodic in the y and z directions. The refined grid covers 
the region x 6 [—0.25, +0.25] . The constraint violation in corresponding unigrid runs (not shown) has the same magnitude. 



E. Gauge Wave 

As another test of the BSSN system we implemented the gauge wave of js^l . From the typical Minkowski coordi- 
nates {t, X, y, z}, one defines new coordinates {t, x, y, z} via the coordinate transformation 

Ad fln{x-t)\ 
Ad (ln{x-t) 



X = X + - — cos ; (19) 

An \ d J ^ ' 

y = y (20) 

z = z, (21) 
where d is the size of the simulation domain. In these new variables, the 4-metric is 

ds^ = -Hdf + Hdx^ + dy^dz^, (22) 

where 

H = H(x-0 = 1-Asin(^^^fc^^. (23) 

This test provides us with an exact solution to which we can compare our numerical results. In addition to the exact 
values a = ^/H and gxx = H, we will compare the extrinsic curvature Kjj, for which the only nonzero component is 

Kxx = --^ COS — ^— '- . (24) 



d^/H 

Since [3' = in the analytic solution we do not evolve the shift but keep it set to zero at all times. 

For this simulation we find it useful to add dissipation to the evolution equations to suppress high frequency 
noise at very high resolutions. The reason is evident from unigrid simulations at high resolutions, as demonstrated 
in Figure IT2I 

The simulation domain was set up in the same way as in Is^l , which is also the same as is used in Section HV DI 
That is, the simulation domain extended almost only in the x direction with x E [—0.5, +0.5]. The domain has only 
3 grid points in the y and z directions. All boundaries are periodic. Although the domain is thus essentially one 
dimensional, the simulation was performed with the full three dimensional code. FieurelT^shows the results after 5 
crossing times for the metric component gxx, comparing refinement and unigrid runs. We see perfect second order 
convergence. 
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FIG. 12: Error in the extrinsic curvature component Kxx along a small part of the x axis for a unigrid simulation after 0.25 crossing 
times. The resolution. Ax = 1/ 2560, is rather high. The noise (that was not present in the initial data) can be significantly reduced 
by dissipation. It is surprising to see that a unigrid simulation with the plain BSSN formulation shows this behaviour; this might 
point to an instability in the system of equations. 
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FIG. 13: Errors for a uniform grid (left) and with mesh refinement (right) after 5 crossing times for the three-metric component 
gxx in the gauge wave test case. The errors have been scaled with the resolution according to second order convergence. 



F. Schwarzschild Black Hole with Excision 



Next, we evolved a static Schwarzschild black hole in Kerr-Schild coordinates which are manifestly static. The line 
element is 

As- = - (l - ^) rf^^ + (^) (l + ^) + (25) 

where M is the mass of the black hole. We choose to use excision to handle the singularity. We use the excision 
method together with one of the tests described in fi^ . The shift is held static at the analytic solution and the lapse 
is evolved using 1 + log slicing. For simplicity we do not search for an apparent horizon but merely excise those 
points within a cube with corners at IM. This small excision region removes the divergences due to the singularity 
while retaining some of the steep gradients. However, as noted by |53, 54], this does not guarantee that the future 
light cone at the excision boundary is contained within the excised region. 

We evolve only one octant of the grid to take advantage of the symmetries present. We use a coarse grid spacing of 
Ax" = 0.4M and At^ = O.IM with 29^ points, giving an outer boundary at lOM (ignoring two symmetry points and 
two outer boundary points). The fine grid also contains 29^ points. The points are not staggered about the origin, so 
there is always a grid point at r = which must be excised. 

We compare the runs using mesh refinement with simulations using the unigrid code as described in f45', 46]. The 
coarse unigrid test is identical to the coarse grid in the simulation using refinement. That is, 29^ points are used 
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with a grid spacing of Ax = 0.4M giving an outer boundary at lOM. For the medium unigrid run the resolution 
is doubled whilst the outer boundary location is held fixed, giving the same effective resolution near the excision 
region as the simulation using refinement. That is, 55^ points are used with a grid spacing of Ax = 0.2M. To compare 
with the high resolution run using refinement we also perform a unigrid run with the same effective resolution. That 
is, 105'^ points are used with a grid spacing of Ax = O.IM. 




1e-20 I 1 1 1 1 

500 1000 1500 2000 

Time (M) 

FIG. 14: The r.m.s. of the change in the lapse with time. As in |45] we see an exponentially damped oscillation as the system 
settles down. However, at sufficiently high resolution an instability sets in. This appears to come from floating point round-off 
error at the initial time and is clearly not caused by mesh refinement. 

Results showing the norm of the change of the lapse in time are shown in Figure Firstly we note that the use 
of refinement combined with an excision boimdary has no qualitative effect on the simulation. As in 14^ the change 
in the lapse shows an exponentially damped oscillation in the absence of instabilities. 

However, in this case we do see an exponentially growing instability which sets in only at the highest resolutions. 
By tracing back the magnitude of this mode we see that it appears to come from floating point round-off error at 
the initial time. It appears in runs both with and without mesh refinement and the growth rate is the same in both 
cases. The origin of this instability is unclear, especially as very similar simulations with the same resolutions were 
shown to be stable in |45]. However, the important point for this paper is that the stability of the simulations are 
independent of the mesh refinement. 

The convergence of the Hamiltonian constraint at f = 500M is shown in figure^] which should be compared with 
Figure 3 in |45]. Second order convergence for the imigrid and mesh refinement simulations are clear away from 
the excision boundary, whilst at the excision boundary the convergence is not so obvious due to the low resolution. 
The error in the mesh refinement runs is comparable to the unigrid runs with the same effective resolution, although 
only near the excision boundary are the full benefits seen. Given the computational resources required, illustrated 
in TableD the benefits of mesh refinement are clear. 



Effective 


Grid size 


Refinement 


Memory 


Real (user) time 


resolution (M) 




levels 


use (MB) 


(s) 


0.4 


25 


1 


58 


18.0 (17.6) 


0.2 


25 


2 


114 


85.5 (84.7) 


0.2 


50 


1 


303 


222.1 (219.4) 


0.1 


50 


2 


695 


957.6 (950.8) 


0.1 


100 


1 


2024 


3337.9 (3309.9) 



TABLE I: Computational resources required to evolve the single black hole with excision to t = 2M with various grid sizes. The 
results using refinement give comparable results to unigrid results with the same effective resolution whilst using approximately 
30% of the resources. These results used a 2.8GHz Intel Xeon processor with the Intel version 7 compilers. 
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FIG. 15: Second order convergence of the Hamiltonian constraint for the Schwarzschild black hole at t = 500M. For the results 
using mesh refinement we plot the results using the finest possible level, causing a "jump" at the interface between the fine 
and coarse grid. We see that the mesh refinement simulations converge to second order except at the excision boundary. The 
simulation with refinement has comparable accuracy as the unigrid run with the same effective resolution near the excision 
boimdary but is less accurate away from this region. 

V. CONCLUSIONS 

We have shown comparisons of results obtained via fixed mesh refinement (FMR) techniques and via standard 
unigrid methods. Ideally, the appropriate use of FMR should produce results very similar in accuracy to the corre- 
sponding unigrid results, while retaining the stability and convergence properties of the unigrid simulations. "Cor- 
responding" in this case implies that the unigrid simulation has the same resolution as the finest grid of the FMR 
simulation. This ideal is only expected to hold for simulations where most of the "interesting behaviour" occurs in 
a limited region which is covered by refined grids, such as in simulations of compact objects. 

We find that particular care must be taken at boundaries between coarse and fine grids, in order to retain stability 
and convergence. In cases where the system contains second spatial derivatives and a multistep time integration 
method is used, we introduce buffer zones to guarantee the convergence of the mesh refinement implementation. 
The need for these buffer zones can be clearly seen in very fine resolution simulations of the wave equation when 
written in "mixed" form (with second spatial derivatives and first derivatives in time). The same effects are seen in 
simulations of gauge waves even at intermediate resolutions when using the BSSN formulation in a similar mixed 
form. 

In order to use parabolic interpolation in time, even from the initial time, we designed a scheme by which data is 
initially evolved forwards and backwards in time. This means that the fine grid evolution can begin with three time 
levels of coarse grid data with which to do the time interpolation. 

We are able to obtain results from FMR simulations which compare favourably with corresponding unigrid sim- 
ulations for the following systems: the wave equation with 1 /r data and excision near the origin, robust stability 
tests, a nonlinear gauge wave, and a Schwarzschild black hole with excision. The robust stability tests indicate that 
numerical noise combined with inter-grid transport operations will not lead to exponential blow ups, and that the 
introduction of refinement boundaries does not produce instabilities. Both the wave equation test and the gauge 
wave test show that proper second order convergence can be obtained irrespective of waves propagating through 
the interfaces between coarse and fine grids. For the case of the Schwarzschild black hole, although we do see an 
instability at high resolutions, this instability is present in the unigrid system as well. At lower resolutions, we see 
stable behaviour showing close agreement between FMR and unigrid results. 

There are few other implementations of mesh refinement in more than one dimension in numerical relativity with 
which to compare. The method of |18] is very close to ours. However, it is an axisymmetric code designed to 
evolve relativistic scalar fields. The code presented in |25|] uses the ADM formalism (which is known to be unsta- 
ble) and exhibits problems at the refinement boundaries. In the implementation of 1 13] which uses Paramesh (|55]), 
all refinement levels are evolved using the same timestep size. This avoids the difficulties that we observe in ob- 
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taining convergence at refinement boundaries. However it may be less efficient, especially where a large range of 
refinement levels is employed. Our implementation provides both the full efficiency of the Berger-Oliger method 
and is a generic interface allowing the simple implementation of other formulations and systems such as relativistic 
hydrod5mamics. 

Finally, we note that the evolutions described in this paper were performed by taking routines for unigrid simu- 
lations of each physical system of interest, and after only slight modifications to these routines, the original unigrid 
driver was exchanged for the "Carpet" FMR driver. The application of FMR techniques to existing imigrid systems 
is thus something which, from a code development viewpoint, can be performed generally. A clear advantage of this 
approach is that most existing analysis tools (such as for wave extraction or apparent horizon finding) will continue 
to work. Thus we invite other researchers to make use of our freely available code to perform their finite-difference- 
based simulations in an FMR setting, and thus achieve higher accuracies with less computing resources than their 
current imigrid simulations may require. 
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APPENDIX A: COUPLING COARSE AND FINE GRID TIME EVOLUTION 



When certain formulations of the equations are evolved with multi-step time integration methods, we find that 

the standard Berger-Oliger approach leads to a loss of convergence at the boundaries of fine grids. We will consider 
only second order methods here. We will argue that the use of Berger-Oliger style prolongation to find boundary 
conditions for refined grids at the intermediate steps causes a second order truncation error after a fixed number of 
time steps. Therefore the global error will be first order for a hyperbolic system where to integrate to a fixed time the 
number of timesteps increases with the grid resolution. 
We first give an example, considering a simple one dimensional example on an infinite domain. The wave equation 

dju = dlu (Al) 

can be reformulated into a system that is first order in time, but still second order in space: 

dtu = V {AT) 
dfV = dlu. (A3) 

Discretising this system on a uniform grid in space and time, we represent the fimction u and v through 

u{t,x) = LT," (A4) 
v{t,x) = Vf (A5) 

at the collocation points t — nAt, x — iAx, where n and i number the grid points in the temporal and spatial direction, 
and At and Ax are the temporal and spatial grid spacings, respectively. The ratio a — At /Ax is the CFL factor, h will 
be used to refer to terms that are order Ax or At. 

For the spatial discretisation of the operator 3^ we will consider the second order centred differencing stencil 



DlF{x) := .^{F{x-Ax)-2F{x)+F{x + Ax)). (A6) 



For time integration we consider one form of the ICN (iterative Crank-Nicholson) method. It consists first of an 
Euler step 

Te[Fo] := Fo + AtdtFo {AT) 
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and then several (in our case, two) ICN iterations 

Tj[Fo,F] := Fo + Atdt^{Fo+F) (A8) 
so that the overall ICN integration with 2 iterations is given by 

T'iCN{2)[fo] := Ti[Fo,Ti[Fo,TEm]. (A9) 

Note that there exist several slightly different variants of the ICN method. 

This discretisation is correct up to second order, which means that for a fixed CFL factor a, the discretisation error 
is 0{h^). This means also that functions f{t,x), which depend on t and x in a polynomial manner with an order 
that is not higher than quadratic, are represented exactly, i.e. without discretisation error. (Such functions / can be 
written as 

fit,x) := t L CpcjtJ'x^ (AlO) 

p=0 q=0 

with constant coefficients Cpq.) 



1. Performing a single step 

It is illustrating to perform a time integration step by step. We will use for that the solution of the wave equation 

u{t,x) = \i^ + \x^ (All) 

v{t,x) = t (A12) 

which will test whether the formulation is indeed capable of representing all functions of the form JAlOt . 
Starting from the exact solution at time to, 

Uito,x) = \il + \x^ (A13) 

y(fo,x) = io, (A14) 
the result of the Euler step of the ICN integration is 

u(°'(fo + Af,x) = U{tQ,x) + MV{tQ,x) (A15) 

= \tl+^-x' + {M)t^ (A16) 

V(°)(f0 + Af,x) = V{tQ,x)+ AtDlU{tQ,x) (A17) 

= fo + Af (A18) 

where F^^^i denotes the result after k ICN iterations. In these expressions, ll^'^\tQ + At,x) differs from the true 
solution U{tQ + h, x) by a term 0{h^). 
The first ICN iteration then leads to 

Jj(i)(fo + At,x) = U{to,x) + At^(v{to,x) + V<^°\to + At,x)^ (A19) 

= ltl + ^x^ + Atto + ^{Atf (A20) 

= l(to + Af)2 + ix2 (A21) 

V^^\to + At,x) = to + At^(DlUito,x) + DlU^°\to + At,x)^ (A22) 

= to + At (A23) 
which is already the correct result. The second ICN iteration does not change the above expressions. 
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2. Mesh refinement 



Let us now introduce mesh refinement, so that the spatial and temporal resolution is not uniform any more. Let 
the grid points be staggered about the origin, so that the grid points with x < have a coarse spatial grid spacing 
lAx, and the grid points with x > have a fine spatial grid spacing Ax. Let the CFL factor a remain uniform, so that 
the temporal grid spacing is correspondingly 2At for x < and Af for x > 0. We will use the Berger-Oliger time 
stepping scheme, which is explained in Section lll Cl and illustrated in FigureOlof the main text. 

The discretisation of 3^. needs so-called ghost points on the fine grid, which are filled by prolongation. Our third 
order polynomial spatial prolongation operator is given by 

{F(x) when aligned 

T^[-l,+9,+9,-l] (A24) 
■F(x+ [-3, -1, +l,+3] Ax) otherwise 

and the second order polynomial temporal prolongation operator by 

{F(t) when aligned 

i[-l,+6,+3] (A25) 
■ F(f+ [-3,-l,+l] At) otherwise 

where [■ ■ ■] denotes a vector, and the operator (■) a dot product. 

Let us now consider a time evolution of the above solution llAllMXT^ . According to the Berger-Oliger time 
stepping scheme, the coarse grid evolution happens as on a uniform grid. The fine grid can now be evolved in two 
different ways: (a) with a Dirichlet boundary condition from prolongating from the coarse grid, or (b) without a 
boimdary condition, i.e. only as an IVP (as opposed to an IBVP). In case (b), the "lost" points have to be filled by 
prolongation after the time step; in that case, prolongation therefore has to fill not 1, but A: + 1 grid points for an ICN 
integration with k iterations. This is illustrated in Figure|4]of the main text. 



3. ICN integration with prolongation 



Let us examine case (a) in more detail, in which case the fine grid boundary values at x < are given through 
prolongation. For the solution ^AlHiAT2t . time integration and prolongation are exact, hence it is for x < 

Lr«(to + Af,x) = l(to + Af)2 + ix2 (A26) 

yW(to + Af,x) = to + At (A27) 

for all ICN iterations fc. (This assumes that all ICN iterations end up the final time. There are also different ways of 
expressing ICN.) 

The Euler step then leads to, for x > 0, 

!j(°)(fo + At,x) = if§ + ix2 + Affo (A28) 
yW(fo + At,x) = h + At (A29) 
as in JAISMaTsI above. For all values of x, this can be written as 

!jW(to + At,x) = i(to + A02 + ix2-0(x>0) i(Ax)2 (A30) 
y(°)(to + At,x) = fo + Af (A31) 
by using the @ fimction, which is defined as 

e(L) = 1^ when Lis false ^^32^ 
1 when L is true. 
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The first ICN iteration leads to 

LZ(i)(io + Ai,x) = l(fo + Af)2 + lx2 (A33) 
y(i'(to + Af,x) = to + At + e{0<x<H) ^a^At. (A34) 

On a uniform grid, and for our solution u, the first ICN iteration is already exact. In general, the first ICN iteration 
should lead to an error which is 0{h?-). This is here not the case; clearly V^^' has an 0{h) error. However, the error 
is confined to a region of length Ax, so an integral norm of the error can be O(ft^); it is not clear what this means in 
practice. 

The result of the second ICN iteration is 

!j(2)(to + Af,x) = i (fo + Af)^ + ^x2 + e(0<x<H) ^^^(Af)^ (A35) 

v'^^\tQ + At,x) = tQ + At. (A36) 

This expression is the final result of the first fine grid time step. It has an error that is 0{h^), or 0{h^) in an integral 
norm. As noted above a local error that is 0{h?-) will lead to a global error of 0{h) when considered at a fixed time. 
That is, the method is only first order convergent. 

A third ICN step — that we don't perform, but it would be possible to do so — would result in 

Lr(3)(fo + At,x) = i(to + Af)2 + ix2 (A37) 

v(3)(fo+2Ai,x) = to + At (A38) 

+ [-2G(0<x<Ax) +0(Ax<x<2Ax)] -a^At 

8 

which is worse than the result of the previous iteration. The error is again 0{h), but has now "infected" two grid 
points. The error is smaller by a factor a^, but the error at the two grid points has a different sign, indicating the start 
of an oscillation. 

It should be noted that using higher order derivatives, e.g. using a five-point stencil for a fourth-order second 
derivative, does not remove this error. Similarly, using a higher order prolongation operator does not help. (We 
assume that it would be possible to adapt the prolongation scheme directly to the time integration scheme and 
arrive at a consistent discretisation in this way. We do not think that this is worthwhile in practice.) 

The main problem seems to be caused by taking a second derivative, which has formally an O ( 1 ) error. The usual 
arguments why the error should be smaller in practice [561 do not hold near the discontinuity that is introduced at 
the refinement boundary. 

The result of simulating case (a) is shown in Figure|7|of the main text. For high resolutions, convergence seems to 
be only of first order. Figure |S1 shows case (b), which shows second order convergence for all resolutions. 



4. The general case 

In the general case the coarse grid evolution and the prolongation will not provide exact data at the boundary on 
the fine grid. Instead we argue heuristically as follows. The initial data U^'^\ V^'^^ will be correct to order h^, 

U" = U{t",x) + }i^e"f"{x) + 0{h^) (A39) 

where e" is a constant and /" a smooth function of x. This data will be evolved forward on the coarse grid using ICN 
to time where it will also be correct to order /z^, 

LI"+^ = iJ(t"+\ x) + /z2e"+i/"+i (x) +0{h^) (A40) 

The fine grid data will then be evolved. The first ICN step is an Euler step to f"+i/2 (as the fine grid timestep is 
one half the coarse grid timestep). On the interior of the refined grid the result will be first order accurate in time 
and second order in space, 

Un+i/i ^ u{t"+y^, x) + {Ate';+^^^ + {Atfe'i+^^^) f+^^x) + 0{h^). (A41) 
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The boundary data for the fine grid is then found by prolongation. Interpolation in space and time, assumed to be 
second order accurate, gives 

!J"+i/2(x < 0) = Lr(f"+^/^ x) + h^el^^'^f"+^l^{x) + 0{h^). (A42) 
Therefore at the refinement boimdary there is a discontinuous jump in the function. 



x=0+ 
x=0- 



he"^^"^ + {Atf (el-^'l^ - e"/'^') ) f"+'l\x). (A43) 



As shown in the case of simple polynomial data in appendix IA 31 the next step in the ICN loop will lead to a second 
order local error when the second derivative of the function is taken at the discontinuity. Thus the error at a fixed 
time will only be first order convergent. It is clear that this discontinuity will also lead to first order errors with other 
multi-step time integration methods such as the Runge-Kutta methods. 

We have performed a symbolic calculation of two complete coarse grid time steps with generic initial data using 
Maple 1 57]. We find that the errors in u and v scale with Ax^ and Ax, respectively, both after the first and the second 
coarse grid steps. This is consistent with the calculation above as well as our numerical results. This means that, 
after a fixed time t, the error in m is 0(Ax), so that we expect first order convergence only. 
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